Mon. Not. R. Astron. Soc. 0Q0Q.[TllT6l(201 1) Printed 27 February 2013 (MN MfiX style file v2.2) 



SubHaloes going Notts: The SubHalo-Finder Comparison Project 

Julian Onions/* Alexander Knebe,^ Frazer R. Pearce/ Stuart I. Muldrew/ 
Hanni Lux/ Steffen R. Knollmann/ Yago Ascasibar/ Peter Behroozi/'^'^ 
Pascal Elahi,^ Jiaxin Han,^'''^ Michal Maciejewski,^ Manuel E. Merchan/'^ 
Mark Neyrinck/^ Andres N. Ruiz/° Mario A. Sgro/'^ Volker Springel/^'^^ 
(N and Dylan Tweed^^ 

T~H ^School of Physics & Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK 

C _ ^ , ^ Departamento de Fi'sica Teorica, Modulo C-15, Facultad de Ciencias, Universidad Autonoina de Madrid, 28049 Cantoblanco, Madrid, Spain 
Cn ' ^ Kavli Institute for Particle Astrophysics and Cosmology, Stanford, CA 94309, USA 
» I ' '^Physics Department, Stanford University, Stanford, CA 94305, USA 
C^ I ^SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA 

^Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Shanghai 200030, China 

^Graduate School of the Chinese Academy of Sciences, 19A, Yuquan Road, Beijing, China 
VQ ' ^Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DHl 3LE, UK 

^ Max-Planck-Institut fur Astrophysik, Garching, Karl-Schwarzschild-Strafie 1, 85741 Garching bei MUnchen, Germany 

^^Instituto de Astronomia Teorica y Experimental (CCT Cordoba, CONICET, UNC), Laprida 922, X5000BGT, Cordoba, Argentina 

^^ Department of Physics and Astronomy, Johns Hopkins University, 3701 San Martin Drive, Baltimore, MD 21218, USA 

^^ Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany 

^^ Zentrum fiir Astronomic der Universitdt Heidelberg, ARl, Monchhofstr. 12-14, 69120 Heidelberg, Germany 

^"^Institut d'Astrophysique Spatiale, CNRS/Universite Paris-Sud 11, 91405 Orsay, France 



6 



27 February 2013 



ABSTRACT 

We present a detailed comparison of the substructure properties of a single Milky Way sized 
dark matter halo from the Aquarius suite (Springel et al.) at five different resolutions, as identi- 
K^ ■ fied by a variety of different (sub-)halo finders for simulations of cosmic structure formation. 

'^ \ These finders span a wide range of techniques and methodologies to extract and quantify 

substructures within a larger non-homogeneous background density (e.g. a host halo). This 
includes real-space, phase-space, velocity-space and time-space based finders, as well as find- 
ers employing a Voronoi tessellation, friends-of-friends techniques, or refined meshes as the 
CO . starting point for locating substructure. A common post-processing pipeline was used to uni- 

^^ ' formly analyse the particle lists provided by each finder We extract quantitative and compara- 

ble measures for the subhaloes, primarily focusing on mass and the peak of the rotation curve 
for this particular study. We find that all of the finders agree extremely well on the presence 
and location of substructure and even for properties relating to the inner part part of the sub- 
halo (e.g. the maximum value of the rotation curve). For properties that rely on particles near 
the outer edge of the subhalo the agreement is at around the 20 per cent level. We find that 
j^ \ basic properties (mass, maximum circular velocity) of a subhalo can be reliably recovered if 

the subhalo contains more than 100 particles although its presence can be reliably inferred for 
a lower particle number limit of 20. We finally note that the logarithmic slope of the subhalo 
cumulative number count is remarkably consistent and < 1 for all the finders that reached 
high resolution. If correct, this would indicate that the larger and more massive, respectively, 
substructures are the most dynamically interesting and that higher levels of the (sub-)subhalo 
hierarchy become progressively less important. 
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1 INTRODUCTION 

The growth of structure via a hierarchical series of mergers is 
now a well established paradigm (_White & Rees 1978). As larger 
structures grow they subsume small infalling objects. However the 
memory of the existence of these substructures is not immediately 
erased, either in the observable Universe (where thousands of indi- 
vidual galaxies within a galaxy cluster are obvious markers of this 
pre-existing structure) or within numerical models, first noted for 
the latter bv lKlvpin etaP ( Il999h . 

Knowing the properties of substructure created in cosmologi- 
cal A'^-body simulations allows the most direct comparison between 
these simulations and observations of the Universe. The fraction of 
material that remains undispersed and so survives as separate struc- 
tures within larger haloes is an important quantity for both studies 
of dark-matter detec tion (Springel et al. 2008: Kuhl en et al.ll200^ : 
IVoeelsberger et al.l2009 : Zavala et al. lOlO") and the apparent over- 
abundance of substru cture within numerical models when com- 
pared to observations dKlvpin et al.ll 19991 : iMoore et alj|l999b . The 
mass and radial position of the most massive Milky Way satellites 
seem to raise new concerns fo r our standard ACDM cosmology 
JBovlan-Kolchin et alJIJOllbl Q: Idi Cintio et ai]|201ll : iFerrero et all 
120111) . while differences between the simulated and observed in- 
ternal density profiles of the satellites seems to have been rec- 
onciled by taking baryonic eff ects into account (e.g. !0h et al.l 
I2OIII : IPontzen & Govematdl201 ih . We are certain that between 5 
per cent and 10 per cent of the material within simula ted galac- 
tic si z ed haloes exists with i n bound substructu res (e.g. iGao et al.l 
I2OO4I : be Lucia et aDl2004l : IContini et ai]l201lh and a substantial 
part of t he host halo has formed from disrupted subhalo mate- 



rial (e.g. iGill et al.l 



ICooper et al.ll2010l : 



2004: Knebe et al. 200^; IWamick et ai]|2008l : 
Libeskind et al.ii201 L) . 



Quantification of the amount of substructure (both observa- 
tionally and in simulations of structure formation) is therefore an 
essential tool to what is nowadays referred t o as "Near-Field Cos- 
mology" uFreeman & Bland-Hawthomll2002h and attempts to do so 
in numerical models have followed two broad approaches: either a 
small number of individual haloes are simulated at exquisite reso- 
lution (e.g. lDiemand et alJl2008l : ISpringel et al.ll2008l : IStadel et al.l 
I2OO9I respectively the Via Lactea, Aquarius and GHalo projects) or 
a larger representative sample of the Universe is mod elled in order 
to quantify halo-to-halo substructure variations (e.g. lAngulo et al.l 
[2009: Bovlan-Kolchin et al. 2009: Klypin et al. 201 1, who used the 
Millennium simulation, Millennium II simulation and the Bolshoi 
simulation, respectively). As this paper studies the convergence of 
halo-finders within a single halo we can add nothing to the topic of 
halo-to-halo substructure variations. 

In a very comprehens ive study that include d 6 different haloes 
and 5 levels of resolution ISpringel et al.l 11200 8h utilised their sub- 
structure finder SUBFIND to detect around 300000 substructures 
within the virial radius of their best resolved halo. They found 
that the number counts of substructures per logarithmic decade 
in mass falls with a power law index of at most 0.93, indicating 
that smaller substructures are progressively less dynamically im- 
portant and that the central regions of the host dark matter halo are 
likely to be dominated by a diffuse dark matter component com- 
posed ofhiuiA«ds_of_thousaiids of streams of tidally stripped ma- 
terial. |Macigiewski_et_alJ(l201jJ) confirmed the existence and prop- 
erties of this stripped material using a 6-dimensional phase-space 
finder HSF. A similar power law index was also found for the larger 
cosm ological studies ( Bovlan-Kolchin et al. 2009: Angulo et al.l 
I2OO9I) . For the Bolshoi simulation iKlvpin et alj J20111) find results 



that are in agreement with their re-analysis of the Via Lactea 11 of 
lOieman d et al. (2008) with an abundance of subhaloes falling as 
the cube of the subhalo rotation velocity. Rather than the present 
value of the maximum rotation velocity they prefer to use the value 
that the subhalo had when it first became a subhalo (i.e. on infall). 
This negates the effects of tidal stripping and harassment within the 
cluster environment but makes it difficult for us to directly compare 
as we have generally only used the final z = Q snapshot for this 
comparison study. 

In recent years there has not only been a number of dif- 
ferent groups performing billion particle single-halo calculations, 
there has also been an explosion in the number of methods avail- 
able for quantifying the size and location of the stru ctures within 
such an A'^-body simulation (see, for instance, Fig.l in lKnebe et alJ 
[201 1). In this paper we extend the halo finder comparison study 
of IKnebe et alj <201ll ) to examine how well these finders extract 
the properties of those haloes that survive the merging process and 
live within larger haloes. Wh ile this issue has already been ad- 
dressed bv iKnebe et alJ <201lh it was nevertheless only in an aca- 
demic way where controlled set-ups of individual subhaloes placed 
into generic host haloes were studied: here we apply the compar- 
ison to a fully self-consistently formed dark matter halo extracted 
from a cosmological simulation. As the results of credible and reli- 
able subhalo identification have such important implications across 
a wide range of astrophysics it is essential to ask how well the 
(sub-)halo finders perform at reliably extracting subhaloes. This 
still leaves open the question of how well different modern grav- 
ity solvers compare when performing the same simulation but at 
least we can hope to ascertain whether or not - given the same set 
of simulation data - the different finders will arrive at the same 
conclusions about the enclosed subhalo properties. We intend this 
paper to form the first part of a series of comparisons. It primar- 
ily focuses on the most relevant subhalo properties, i.e. location, 
mass spectrum and the distribution of the value of the peak of the 
subhaloes' rotation curve. 

In Section |2] we begin by summarising the eleven substruc- 
ture finders that have participated in this study, focusing upon any 
elements that are of particular relevance for substructure finding. 
In Section [3] we introduce the Aquarius dataset that the described 
finders analysed for this study. Both a qualitative and a quantitative 
comparison between the finders is contained in Section |4] which 
also contains a discussion of our results, before we summarise and 
conclude in Section [S] 



2 THE SUBHALO FINDERS 

In this section we present the (sub-)halo finders participating in the 
comparison project in alphabetical order. Please note that we pri- 
marily only provide references to the actual code description pa- 
pers and not an exhaustive portrait of each finder as this would be 
far beyond the scope of this paper. While the general mode of op- 
eration can be found elsewhere, we nevertheless focus here on the 
way each code collects and defines the set of particles belonging 
to a subhalo: as already mentioned before, those particle lists are 
subjected to a common post-processing pipeline and hence the re- 
trieval of this list is the only relevant piece of information as far as 
the comparison in this particular paper is concerned. 
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Figure 1. The images show the smoothed dark matter density within 
of the recovered subhaloes for the finder labelled at the top of each 
than 10 km/s are shown. 
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dark matter density within an quadrant at resolution level 4. In each panel the overplotted circles indicate the location 
labelled at the top of each panel. They are scaled proportionally using Umax ■ Only subhaloes with a iimax greater 
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Figure 1. (continued) Recovered subhalo locations and iimax scale by labelled finder. 
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2.1 ADAPTAHOP (Tweed) 

ADAPTAHOP is a full topological algorithm. The first stage consist 
in estimating a local density using a 20 particles SPH kernel. Par- 
ticles are then sorted into groups around a local density maximum. 
And saddle points act as links between group. All groups are first 
supposed to be one single entity, that we hierarchically divide into 
smaller groups, by using an increasing density threshold. Haloes 
are then defined as groups of groups linked by saddle points cor- 
responding to densities higher that 80 times the mean DM density. 
By increasing the threshold we further detail the structure of the 
halo as a node structure tree. Where a node is either a local max- 
ima, or a group of particules connecting higher level nodes. After 
using this bottom to top approach, the (sub)haloes are defined us- 
ing a top to bottom approach, hierarchically regrouping nodes so 
that a sub(sub)halo has a smaller mass than its host (sub)halo. Each 
particle belongs to a single structure either a halo or a subhalo. The 
(sub)haloes' centers are defined as the position of its particles with 
the highest SPH density. We need to stress that no unbinding pro- 
cedures are used in this algorithm, at the risk of overestimating the 
number/misidentification of subhaloes with a low number of parti- 
cles. 



2.2 AHF (Knollmann & Knebe) 

The halo finder AWv^ (AMIGA Halo Finder) is a spherical over- 
d ensity finder that identifies (isolated and sub-)haloes as described 
in lGilletalJj2004h as well as lKnoUmann & Knebd t200gl) . The ini- 
tial particle lists are obtained by a rather elaborate scheme: for each 
subhalo the distance to its nearest more massive (sub-)halo is cal- 
culated and all particles within a sphere of radius half this distance 
are considered prospective subhalo constituents. This list is then 
pruned by an iterative unbinding procedure using the (fixed) sub- 
halo centre as given by the local density peak determined from an 
adaptive mesh refinement hierarchy. For more details we refer the 
reader to aforementioned code description papers as well as the on- 
line documentation. 



2.3 Hierarchical Bound-Tracing (HBT) (Han) 

HBT JHan et al J 1 20 111) is a tracing algorithm working in the time 
domain of each subhaloes' evolution. Haloes are identified with 
a Friends-of-Friends algorithm and halo merger trees are con- 
structed. HBT then traverses the halo merger trees from the earli- 
est to the latest time and identifies a self-bound remnant for every 
halo at every snapshot after infall. Care has been taken to ensure 
that subhaloes are robustly traced over long periods. The merg- 
ing hierarchy of progenitor haloes are recorded to efficiently allow 
satellite-satellite mergers or satellite accretion|j 



2.4 HOT-i-riEstAS (HOT3d & hot6d) (Ascaslbar) 

HOT-l-FiEstAS is a general-purpose clustering analysis tool, still 
under development, that performs the unsupervised classifica- 
tion of a multidimensional data set by computing its Hierarchi- 
cal Overdensity Tree (HOT), analogous to the Minimal Span- 
ning Tree (MST) in Euclidean spaces, based on the density 



^ AHF is freely available from http://www.popia.ft.uam.es/AMIGA 

^ It should be noted that HBT had access to the full snapshot data for 

Aquarius-A. 



field returned by the Field Estimator for Arbitrary Spaces (FiEs- 
tAS Ascasibar & Binn ev 2005; Ascasibar 2010). As explained in 
lKnebeetalJi201ll) in the context of halo finding, HOT-l-FiEstAS 
identifies objects with density maxima, either in configuration 
space (considering particle positions alone, H0T3d) or in the full, 
six-dimensional phase-space of particle positions and velocities 
(HOT6d). In both cases, the boundary of an object is always set 
by the isodensity contour crossing a saddle point, and its centre is 
defined as the density-weighted average of its constituent particles. 
The main d ifference with respect to the version used in 
I Knebe et alj (2011 ) is that there is now a post-processing stage, akin 
to a 'hard' expectation-maximization that is specifically tailored to 
the problem of halo finding, where: 

(i) ''max and ?;max are computed for every object in the catalog. 

(ii) Objects with more than 10 particles within rmax are labelled 
as (sub)-halo candidates. 

(iii) Particles are assigned to the candidate that contributes most 
to the phase-s pace density at th eir location, approximating each 
candidate by a iHemguistl ( 1 19901) sphere with the appropriate val- 



ues of rn 



and Un 



Candidates are only kept if they contain more than 5 particles 
within rmax and the density within that radius is higher than 100 
times the critical density. Although a detailed discussion is obvi- 
ously beyond the scope of this work, it is interesting to comment 
that some of the objects discarded by the latter criterion seem to be 
numerical artefacts, but others are clearly associated to filaments, 
streams, and other loose - yet physical, sometimes even gravita- 
tionally bound - structures. Since they are certainly not individ- 
ual dark matter (sub)-haloes, they can be simply discarded for our 
present purposes. 



2.5 Hierarchical Structure Finder (HSF) (Maciejewskl) 

The Hierarchical Structure Finder (HSF) identifies objects as con- 
nected self-bound particle sets above some density threshold. This 
method has two steps. Each particle is first linked to a local dark 
matter phase-space density maximum by following the gradient 
of a particle-based estimate of the underlying dark matter phase- 
space density field. The particle set attached to a given maximum 
defines a candidate structure. In a second step, particles which 
are gravitationally unbound to the structure are discarded until 
a fully self-bound final object is obtained. For more details see 
iMacieiewski et alJ ( l2009l) . 



2.6 MENDIETA (Sgro, Ruiz & Merchan) 

The MENDIETA finder is a Friends-of-Friends based finder that is 
used to obtain a dark matter halo. This prospective host halo is 
subsequently refined by looking at peaks of increasing density by 
reducing the linking length. This approach decomposes the halo 
into its substructure plus other minor overdensities. In a final pass 
pass unbound particles are removed by check ing their associat ed 
energies. MENDIETA is described more fully in lSgro et aTl fcOlOl) . 



2.7 ROCKSTAR (Behroozi) 

ROCKSTAR (Robust Overdensity Calculation using K-Space 
Topologically Adaptive Refinement) is a phase-space halo 
finder designed to m aximize halo consistency across timesteps 
JBehroozi et al.ll201 m . The algorithm first selects particle groups 
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with a 3D Friends-of-Friends variant witli a very large linking 
length (6 = 0.28). For each main FOF group, Rockstar builds a 
hierarchy of FOF subgroups in phase-space by progressively and 
adaptively reducing the linking length, so that a tunable fraction 
(70%, for this analysis) of particles are captured at each subgroup 
as compared to the immediate parent group. When this is complete, 
Rockstar converts FOF subgroups into seed haloes beginning at the 
deepest level of the hierarchy. If a particular group has multiple 
subgroups, then particles are assigned to the subgroups' seed haloes 
based on their phase-space proximity. This process is repeated at all 
levels of the hierarchy until all particles in the base FOF group have 
been assigned to haloes. Unbinding is performed using the full par- 
ticle potentials; halo centres and velocities are calculated in a small 
region close to the phase-space density maximum. 



2.8 STF (Elahi) 

The STructure Finder Hierarchical Structure Finder (Elahi et al.l 
I2OIII (STF)) identifies objects by utilizing the fact that dynami- 
cally distinct substructures in a halo will have a local velocity dis- 
tribution that differs significantly from the mean, i.e. smooth back- 
ground halo. This method consists of two main steps, identifying 
particles that appear dynamically distinct and linking this outlier 
population using a Friends-of-Friends-like approach. Since this ap- 
proach is capable of not only finding subhaloes, but tidal streams 
surrounding subhaloes as well as tidal streams from completely dis- 
rupted subhaloes, we also ensure that a group is self-bound. Parti- 
cles which are gravitationally unbound to a candidate subhalo are 
discarded until a fully self-bound is obtained or the object consists 
of fewer than 20 particles, at which point the group is removed en- 
tirely. 



2.9 Subflnd (Springel) 

SUBFIND identifies substructures as locally overdense, gravitation- 
ally bound groups of particles. Starting with a halo identified 
through the Friends-of-Friends algorithm, a local density is esti- 
mated for each particle with adaptive kernel estimation using a pre- 
scribed number of smoothing neighbours. Starting from isolated 
density peaks, additional particles are added in sequence of de- 
creasing density. Whenever a saddle point in the global density 
field is reached that connects two disjoint overdense regions, the 
smaller structure is treated as a substructure candidate, followed 
by merging the two regions. All substructure candidates are sub- 
jected to an iterative unbinding procedure with a tree-based calcu- 
la tion of the potential. T he SUBFIND algorithm is discussed in full 
in lSpringeletaDbOOlll . 



2.10 VOBOZ (Neyrinck) 

VOBOZ (INevrinck et al.ll2005h was developed to have little depen- 
dence on free parameters. Density peaks are found using a Voronoi 
tessellation, which gives an adaptive, parameter-free estimate of 
each particle's density and set of neighbours. Each particle is joined 
to the peak that lies up the steepest density gradient from that parti- 
cle. A halo associated with a high density peak (which is defined as 
the VOBOZ centre of the halo) will typically contain smaller density 
peaks. The significance of a halo is judged according to the ratio of 
its central density to a saddle point joining the halo to a halo with 
a higher central density, comparing to a Poisson point process. For 
this project, we impose a 4-cr significance threshold on subhaloes. 



Particles not gravitationally bound to each halo are iteratively re- 
moved, by comparing their potential energies (measured as sums 
over all other particles) to kinetic energies with respect to the veloc- 
ity centroid of the halo's core (i.e. the particles that directly jump up 
density gradients to the peak). In the unbinding process, the least- 
bound particles are removed first; for each halo, the boundedness 
threshold reduces by a factor of \/2 at each iteration, until it reaches 
its true value. 



3 THE DATA 

The data used for this paper forms part of the Aquarius project 
JSprineel et al.ll2008) . It consists of multiple dark matter only re- 
simulations of a Milky Way like halo at a variet y of resolution s 
performed using GADGETS (based on GADGET2.'Springell l2005h . 
We have used the Aquarius-A halo dataset at 2 = for this project. 
This provides 5 levels of resolution, varying in complexity from 
the 2.3 million particles of the lowest resolution (i.e. level 5), up 
to the 4.25 billion particles of the highest resolution (i.e. level 1), 
as shown in Table [T] The underlying cosmology for the Aquarius 
simulations is the same as that used for the Millennium simula- 
tion JSprineel etal.ll2005h i.e. Om = 0.25, t^A = 0.75, erg = 
0.9, ris = l,/i = 0.73. These parame ters are consistent with the 
latest WMAP data JJarosik et aDl201 ih although erg is a little high. 
All the simulations were started at an initial redshift of 127. Precise 
d etails on the set-up and performance of these models can be found 
in lSpringele"taD(l2008h . 

The participants were asked to run their subhalo finders on 
the supplied data and to return a catalogue listing the substructures 
they found. Specifically they were asked to return a list of uniquely 
identified substructures together with a list of all particles associ- 
ated with each subhalo. 

Finders were initially run on the smallest dataset, the Aq-A- 

5 data. This allowed for debugging of the common output foiTnat 
required by the project and some basic checks on the internal con- 
sistency of the data returned from each participant. Once this had 
been achieved each participant scaled up to the higher resolution 
datasets, continuing until they reached the limits of their finder 
and/or the computing resources readily available to them. A sum- 
mary of the number of subhaloes found by each subhalo finder at 
the various levels is contained in Table |2] as well as the size of the 
largest subhalo at level 4. All of the finders that participated in this 
study completed the analysis of the level 4 dataset which is used 
for the main comparison that follows and contains around 6 million 
particles within the region considered, a sphere of radius 250 kpc//i 
around a fiducial centreQ Three of the finders (AHF, ROCKSTAR 

6 SUBFIND) completed the analysis of the very computationally 
demanding level 1 dataset. In addition to these HBT and HSF com- 
pleted level 2 which contains around 160 million particles within 
the region examined here. 

Both the halo finder catalogues (alongside the parti- 
cle ID lists) and our post-processing software (to be de- 
tailed below) are publically available from the web site 
,http://popia.ft.uam.es/SubhaloesGoingNotts^under the Tab "Data". 



^ We adopted a fixed and unique position for the host halo of x = 
57060.4, y = 52618.6, z = 48704.8 kpc/h independent of the resolu- 
tion. 
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Figure 2. Subhalo recoveiy as a function of resolution. Location and size of recovered substructure from level 3 to level 1 for the three finders that reached this 
level. In all panels subhaloes with t)inax> 10 km/s are shown, scaled by Vmax as in Figure^and the background image is the smoothed dark matter density 
at that level. The relevant finder and level are labelled in the top right of each panel. The biggest change between levels is the additional small scale power 
moving the substructure locations. 



4 THE COMPARISON 

We are going to primarily focus on comparing the location of sub- 
haloes (both visually and quantitatively), the mass spectrum, and 
the distribution of the peak value of the rotation curve. The compar- 
ison, however, is based solely upon the provided particle lists and 
not the subhalo catalogues as the latter are based upon each code's 
own definitions and means to determine aforementioned properties 
and hence possibl y introducing "noise" into the comparison (cf. 
iKnebe et alJl201 ih . In order to achieve a fair comparison between 
the respective finders we produced a single analysis pipeline which 
we used to post-process the particle lists provided by each partici- 



pating group. This ensured consistency across our sample while re- 
moving differences due to the adoption of different post-processing 
methodologies and the particular choice of threshold criteria. The 
comparison detailed in this paper is restricted to this uniform post- 
processed dataset. We intend to explore differences due to different 
methodologies in a subsequent work. However, we stress at the out- 
set that our particular chosen post-processing methodology is not 
intended to be unique nor do we put it forward as the best way of 
defining a subhalo. Rather we use a single methodology so that we 
can first answer the most fundamental question: if we agree on a 
single subhalo definition do the different finders agree on the most 
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Table 1. Summary of key numbers for each Aquarius level, the dataset used for this study. A^^ is the number of particles with the highest resolution (lowest 
individual mass). Ni the number of low resolution particles - the sum of the remainder. N2h0 is the number of high resolution particles found within a sphere 
of radius 250 kpc/h from the fiducial centre at each resolution (i.e. those of interest for this study). Mp is the mass of one of these pailicles (in Mq/H). S is 
the resolution increase (mass decrease) for each level relative to level 5. and Sp is the resolution increase relative to the previous level. All particles are dark 
matter particles. 
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Table 2. The number of subhaloes containing 20 or more particles and centres within a sphere of radius 250kpc/h from the fiducial centre found by each finder 
after standardised post-processing (see Section l4~n . Three finders (AHF, ROCKSTAR& SUBFIND) returned results from the highest resolution (level 1) within 
the timescale of this project. Below this we list the number of particles contained within the largest subhalo after post-processing. 



Number of subhaloes within 250kpc/h of the fiducial centre after post processing. 
AHF HBT H0T3D H0T6D HSF MENDIETA ROCKSTAR STF SUBFIND VOBOZ 



Name 



ADAPTAHOP 
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353 
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58 
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2497 
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1493 
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Aq-A-3 


- 


11213 


11693 


- 


- 


11240 


10948 


11797 


10250 


10094 


13343 


Aq-A-2 


- 


38441 


39703 


- 


- 


35445 


- 


38489 


- 


33135 


- 


Aq-A-1 


- 


226802 


- 


- 


- 


- 


- 


235819 


- 


221229 


- 



Number of particles in the largest subhalo within 250kpc/h of the fiducial centre after post processing. 
49076 77225 66470 69307 61581 73167 48387 78565 56990 50114 54685 



Aq-A-4 



fundamental properties they recover? Perhaps surprisingly we will 
see that the answer to this question is broadly yes. 

We did not consider in this paper efficiency of processing, as 
to make a fair comparison the codes would need to run on compara- 
ble machines with a set amount of memory and processors. In this 
instance the finders were run with the resources that were available 
to each of the participants. Some indication of the capabilities of 
the respective finders may be deduced from Table |2] 



4.1 Post-processing pipeline 

Some finders (e.g. AHF) include the mass (and particles) of a sub- 
halo within the encompassing host halo whereas others do not (e.g. 
SUBFIND), preferring each particle to only be associated with a sin- 
gle structure. Either of these approaches has its pros and cons. For 
instance, keeping the subhalo mass as part of the halo mass makes 
it straightforward to calculate the enclosed dynamical mass of any 
object. However, such an approach easily leads to multiple counting 
of mass, particularly if there are many layers of the substructure hi- 
erarchy. In principle though it is not difficult to transform from one 
definition to the other given knowledge of both the halo and particle 
locations. In our study, 5 of the 1 1 finders chose to include the mass 
of subhaloes whereas the other 6 did not. Following our principle 
of creating a uniform analysis pipeline we processed all the parti- 
cle lists to ensure that a particle could only reside within a single 
structure. To this end, we first sorted the returned halo catalogue 
into mass order. Then starting from the smallest halo we performed 
the centring, trimming and overdensity checks detailed below to 
trim the subhalo uniformly. We then tagged the particles contained 
within this object as being within a subhalo before continuing to 
the next largest subhalo and repeating the procedure ignoring par- 
ticles already tagged as being used before. This preserved the max- 



imum depth of the subhalo hierarchy while ensuring that a particle 
could only reside within a single subhalo. We should remark that in 
practice excising all the sub-subhaloes from each subhalo's particle 
list made little difference to any of the results presented here as at 
any level of the subhalo hierarchy only around 10 per cent of the 
material is within a subhalo of the current halo. So sub-subhaloes 
contribute only around 1 per cent of the halo mass, although it can 
affect other properties such as the centre of mass. 

All the particles belonging to the list each finder identified as 
being associated with a subhalo were extracted from the original 
simulation data files to retrieve each particle's position, velocity 
and mass. From this data the centre of mass was first calculated, 
before being refined based on consideration of only the innermost 
10 per cent of these particles, sorted with respect to the initial cen- 
tre of mass. This procedure was repeated until a stable centre was 
found, i.e. until the change in the position was below the actual 
force resolution of the simulation. Once the centre had been defined 
the particles were ordered radially from this point and a rotation 
curve GM{< r)/r and overdensity M{< r)/(47rr'^/3) calculated 
until it dropped below 200 times the critical density pcrit defining 
the subhalo radius R200 and mass M200 ■ AH particles outside -R200 
were removed which was essential in particular for the phase-space 
finders who also considered already stripped material as still being 
part of and belonging to the subhalo. Please note that our post- 
processing pipeline does explicitly not feature an unbinding proce- 
dure as this already formed part of most halo finding algorithms. 
At this point the maximum circular velocity, Umax was obtained by 
smoothing the rotation curve and locating its maximum by search- 
ing both inwards and outwards for a peak in the rotation curve and 
taking the average of these two measures, a process that stabilises 
the measure if the rotation curve is very flat or noisy. 

We emphasise that the precise subhalo properties are some- 
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what sensitive to the definition of the halo centre. Various groups 
use the centre-of-mass as the centre of all material enclosed within 
the subhalo's radius (both with and without including substructure), 
the centre-of-mass of some smaller subset (as here for example), 
the location of the most bound particle, the location of the densest 
particle or the minimum of the gravitational potential. Additionally 
different groups use different methodologies for deciding whether 
or not a particle is bound to a halo as this involves some decisions 
about the global potential and can be a very time consuming pro- 
cess if done fully generally and iteratively. 

Finally the choice of where to place the subhalo edge is also 
problematic. By definition the subhalo resides within some in- 
homogeneous background density and so at some point particles 
cease to belong to it and should rather be associated with the back- 
ground object. Different groups split the host halo from the subhalo 
in different ways and there is no correct method. Without a uniform 
choice these differences can swamp any differences due to actually 
finding subhaloes or not. We stress that our post-processing (where 
we treat each subhalo in isolation) can only remove particles from 
the original list of those particles associated with a subhalo. We 
have therefore tested whether or not our results are sensitive to our 
choice of 200 as an overdensity parameter by re-running our anal- 
ysis with a tighter threshold of 500. Other than making all the sub- 
halo masses smaller this has no noticeable effect on the scatter of 
the cumulative number counts. We therefore decided to stick to the 
original choice of -R200 and M200, respectively. Further, throughout 
the subsequent comparison only haloes with more than 20 (bound) 
particles within i?2oowere used, although some finders detected and 
returned haloes with less particles. 

To summarise, our uniform post-processing pipeline involved 
the following steps, applied iteratively where necessary: 

• The subhalo catalogues were sorted into mass order. 

• Starting from the smallest subhalo, the particles associated 
with the cuiTent subhalo were obtained from the simulation data. 

• Only particles tagged as "not used before" were considered. 

• The centre-of-mass was iteratively calculated using the inner- 
most 50 per cent of particles. (Originally we used the innermost 10 
per cent but found some of the more dispersed sub structures did 
not converge with this value). 

• A value for -R200 was calculated based on an enclosed over- 
density of 200 times the critical density. 

• The subhalo mass and rotation curve peak «max were com- 
puted based on particles inside -R200 • 

• Only substructures containing more than 20 particles were re- 
tained. 



4.2 Visual comparison 

A visual representation of the location and size (based on t^max) 
of the recovered subhaloes at Aquarius level 4 from each of the 
finders is shown in Figure [T] A smoothed colour image of the un- 
derlying dark matter density based on all particles from the original 
Aquarius data are shown in one quadrant of the main halo, and this 
is overplotted with the recovered subhaloes from each finder indi- 
cated by circles whose size is scaled according to Umax (specifically 
Umax (in km/s) divided by 3). This allows a visual comparison be- 
tween the finders. Only haloes with Wmax > 10/i~^km/s are shown. 
We immediately see that most of the finders are very capable of 
extracting the locations of the obvious overdensities in the under- 
lying dark matter field. Wherever you would expect to find a sub- 
halo (given the background density map) one is indeed recovered. 



This demonstrates that substructure finders should be expected to 
work well, recovering the vast majority of the substructure visible 
to the eye. Additionally, if our aforementioned post-processing is 
applied the quantitative agreement between the finders is also ex- 
cellent, with the extracted structures having very similar properties 
between finders (see below). The majority of the finders agree very 
well, reliably and consistently recovering nearly all the subhaloes 
with maximum circular velocities above our threshold. 

While Figure [T] illustrates the agreement between the finders 
at a single Aquarius level (in this case level 4, which all the par- 
ticipating finders have completed), in Figure|2]we construct a sim- 
ilar Figure to illustrate the agreement between levels. We show the 
same quadrant at level 3 to level 1 for the three finders that have 
completed the level 1 analysis (i.e., AHF, ROCKSTAR & SUBFIND); 
we deliberately omitted level 5 and level 4 as the former is not very 
informative and the latter has already been presented in Figure [T] 
As can be seen, the main difference between the different levels 
is in the exact location of the substructures. This changes because 
additional power was added to the Aquarius initial power spectrum 
to produce the additional small objects that foiTn as the resolution 
is increased (fundamentally, the Nyquist frequency has changed as 
there are more available tracers within the higher resolution box). 
This extra power moves the substructure around slightly, and these 
differences are amplified in the, by definition, non-linear region of 
a collapsed object. Despite this the ready agreement between the 
three finders at any single level is clear to see and this is similarly 
true for both the other finders (HBT, HSF) that completed level 2. 
We do not explore the effect of changing the resolution on sub- 
halo extraction in more detail here because that is not the main 
point of this paper, which focuses on how well different finders ex- 
tract substructure relative to each other. Also, this topic has already 
been well studied for SUBFIND using this same suite of models by 
JSprineel et al.l2008h . 



4.3 Subhalo Mass Function 

4.3.1 Level 4 

Perhaps the most straightforward quantitative comparison is sim- 
ply to count the number of subhaloes found above any given mass. 
For Aquarius level 4 this produces the cumulative mass plot (based 
on Af20o) shown in Figure|3] Results from each participating finder 
are shown as a line of the indicated colour. Generally the agree- 
ment is good, with some intrinsic scatter and a couple of outliers 
(particularly ADAPTAHOP and MENDIETA) which do not appear to 
be working as well as the others, finding systematically too many 
or too few subhaloes of any given mass respectively. For ADAPTA- 
HOP we like to remind the reader that this code does not feature a 
procedure where gravitationally unbound particles are removed; we 
therefore expect lower mass haloes stemming from Poisson noise in 
the background host halo to end up in the halo catalogue as well as 
haloes to have a higher mass in general possibly explaining the dis- 
tinct behaviour of this code. But typically the scatter between codes 
is around the 10 per cent level except at the high mass end where it 
is larger as each finder systematically recovers larger or smaller 
masses in general. We like to remind the reader again that this 
scatter is neither due to the inclusion/exclusion of sub-subhaloes 
(which has been taken care of by our post-processing pipeline) nor 
the definition of the halo edge: as the 10 per cent differences still 
remain if choosing -R500 as the subhalo edge. 

Table[2]lists the number of subhaloes found that contain 20 or 
more particles after the uniform post-processing procedure detailed 
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Figure 3. Cumulative number count of subhaloes above the indicated mass found (M200) within a radius of 250 kpc/h from the fiducial halo centre after 
standardised post-processing at resolution level 4 (see Section|4T|for details). The bottom plot shows the relative offset from the mean of the cumulative mass 



above had been performed and within 250 kpc/h of the fiducial cen- 
tre of the main Aquarius halo at each level completed for all the 
eleven finders that participated. These number counts are generally 
remarkably consistent, again with a few outliers as expected from 
Figure |3] The majority of the finders are recovering the substruc- 
tures remarkably well and consistent, respectively. 

As an additional quantitative comparison we list the number 
of particles associated with the largest substructure found by each 
of the finders as the last row of Table |2] All the finders recover a 



structure containing 60, 000 particles ±20 per cent. As Figure |3] 
has shown there is a lot of residual scatter for the highest mass 
haloes even when a uniform post-processing pipeline is used. This 
is most likely due to the different unbinding algorithms used in the 
initial creation of the substructure membership lists which are par- 
ticularly uncertain for these large structures. At the other end of 
the substructure mass scale we have chosen to truncate our com- 
parison at subhaloes co ntaining 20 particle s as this was shown to 
be the practical limit in lKnebe et alj <201lh . Some participants re- 



© 2011 RAS, MNRAS QOQO.fTlfTel 



SubHalo-Finder Comparison 1 1 



10 



11 



n = 0.90 
n = 0.97 
n = 0.91 
n = 0.93 
n = 0.90 







A 



^ 



_AHF (A-5): ao=2.02e + 09; 
AHF (A-4): ao=9.27e + 09; 
.AHF {A-3): ao=3.08e + 09; 
.AHF (A-2): ao=4.B2e + 09; 
.AHF (A-1): ao=2.00e + 09; 

_ Rockstar (A-5): ao=1.33e+10; n = 0.99 
Rockstar (A-4): ao=7.98e + 09; n = 0.96 
. Rockstar (A-3): ao=2.77e + 09; n = 0.90 
.Rockstar (A-2): ao=4.19e + 09; n = 0.92 
.Rockstar (A-1): ao=3.43e + 09; n = 0.89 

SUBFIND (A-5): ao=2.18e + 09; 

SUBFIND (A-4): ao=1.19e+10; 

SUBFIND (A-3): ao=5.10e + 09; 
.SUBFIND (A-2): ao=4.61e + 09; 
.SUBFIND (A-1): ao=1.74e + 09; 




subhaloes inside 250.0 kpc/h-sphere about fiducial centre 
(ao and n from fit to N{>M)=aQM~^) 



10" 



10° 



10'' 



10' 



10" 



10=^ 



10* 



A^aoo [Mg/^] 



Figure 4. Cumulative subhalo mass function (multiplied by M to compress the vertical dynamical range) for all five Aquarius levels for the AHF, ROCKSTAR, 
and SUBFIND finder. We fit the function A'"(> M)/Ntot = ao X M~" between the mass equivalent to 100 particles at each level and lO^MQ/h. Note: the 
data has not been shifted for clarity but is as plotted. 



turned haloes smaller than this as this is their normal practice. They 
all stress that such small subhaloes should be treated with extreme 
caution but that there does appear to be a bound object at these lo- 
cations even if its size is uncertain. We have removed them here for 
the purposes of a fair comparison. 

Other plots that could be considered are those comparing the 
number of subhaloes against radial distance, or fractional mass 
against radial distance. Both these were produced and considered, 
but did not give any further insight to the comparison. 

4.3.2 All Levels 

Cumulative subhalo number counts like that shown for level 4 in 
Figure |3] can be calculated for all completed levels and compared. 
As shown in Figure |2] while increasing the resolution does not ex- 
actly reproduce the same substructures a reasonable approximation 
is achieved and so we expect to find a set of similar subhaloes con- 
taining more particles as we decrease the individual particle mass 
between levels (i.e. any specific subhalo should effectively be better 
resolved as the resolution increases). We show the cumulative num- 



ber counts for the finders AHF, ROCKSTAR, and SUBFIND (multi- 
plied by M to compensate for the large vertical scale) from level 5 
to level 1 in Figure |4] We show this as an example and stress that 
similar plots with similar features could be produced for any of the 
finders that completed level 2. The curve for each level starts at 20 
particles per halo and we like to stress that no artificial shifting has 
been applied: any differences seen in the plot are due to the dif- 
ferent halo finding algorithms. Below about 100 particles per halo 
the cumulative number counts fall below the better resolved curves, 
indicating that subhaloes containing between 20 and 100 particles 
are not fully resolved and should have a sli ghtly higher associated 
mass, also reported in lMuldrew et al.l ( l201ll) . Above lO^MQ/h the 
power law slope breaks as there are less than 10 subhaloes more 
massive than this limit and the number of these is a property of this 
particular host halo. For these reasons we fit a power law of the 
form 



N{> M) 
Ntot 



■ aoM 



(1) 
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between 100 particles and IO^Mq/Ii where the power law breaks. 
Here ao is a normalisation (capturing the rise in the number of sub- 
haloes due to the increase in resolution), M is the mass and n is the 
power law slope. The fitted values of the parameters by level are 
given in the legend for each finder. The subhalo cumulative num- 
ber count appears to be an unbroken power law - at least in the 
range co nsidered for the fittin g. Similar results for SUB FIND were 
found bv lSpringel etal] ( l2008h . 

We extended this particular analysis of fitting a single power- 
law to the (cumulative) subhalo mass function to all finders at all 
available levels and compare the values of ao and n as a function 
of level for all participating substructure finders in Figure |5] There 
we find that at level 5 little can be said because the fitting range 
is so narrow. At the lower, better resolved levels good agreement is 
seen between the finders (clearly ADAPTAHOP is a strong outlier on 
this plot, probably due to its lack of unbinding as mentioned before 
when discussing Figure[3] and HOT3d (as well as the first resolu- 
tion step of MENDIETA) is inverted with respect to the main trend) 
and a consistent trend emerges: all agree that the power law slope, n 
is less than 1 and if anything decreasing with increasing resolution. 
Values of n less than 1 are significant because they imply that not 
all the mass is contained within substructures, with some material 
being part of the background halo. This has important ramifica- 
tions for studies requiring the fraction of material within substruc- 
tures such as the dark matter annihilation signal and lensing work. 
Although this result is robust between all high-resolution finders 
we remind the reader that this is for a single halo within a sin- 
gle cosmological model. However, it does indicate that, as perhaps 
expected, the most important contribution to substructure mass is 
from the most massive objects and that progressively smaller struc- 
tures contribute less and less to the signal. 



4.4 Distribution of Umax 

If, instead of quantifying the total mass of each subhalo, we rather 
use the maximum rotational velocity, iimax to rank order the sub- 
haloe s in size we obtain a generally much tighter relation (see be- 
low). [Knebe_et_alJ (20n|) already found that Umax was a particu- 
larly good m etric for comparing h aloes and we confirm this for 
subhaloes. As lMuldrew et ahUlOllh show ed in figure 6 of thei r pa- 
per, this is because for an NEW profile (INavarro et al.lll997l) the 
maximum of the rotation curve is reached at less than 20 per cent 
of the virial radius for objects in this mass range so iimax is a prop- 
erty that depends upon only the very inner part of the subhalo and is 
not affected by any assumptions mad e about the outer edge. On the 
other hand, it has also been shown (lAscasibar & Gottlobeill200a) 
that iimax provides a meaningful tracer of the depth of the gravita- 
tional potential (i.e. the mass scale) of the halo. 

Figure |6] displays the cumulative Umax for all the finders for 
level 4 again. All the finders align incredibly well for the largest 
subhaloes with Umax> 20 km/s. For subhaloes smaller than this the 
alignment remains tighter than the total mass comparison down to 
rotation velocities of around 6 km/s. At level 4 haloes of this size 
contain around 80 particles in total, so «max is being calculated 
from less than 20 particles at this point; the arrows give an indica- 
tion of the number of particles inside rmax. ADAPTAHOP, despite 
its missing unbinding procedure, agrees well with other finders for 
high rotation velocities as this particular statistic probes inner re- 
gions of the subhaloes which are less affected by unbound particles; 
and its deviation at the lower Umax -end is due to the existence of 
(small mass) fluke objects not removed by such an unbinding step. 
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Figure 5. A comparison of the slope and normalisation of the fits of the 
mass function derived as per Figure|4]for all finders at all levels returned. 



4.5 Radial Mass Distribution 

The accumulated total mass of material with subhaloes is measured 
by ordering the subhalo centres in radial distance from the fiducial 
centre of the halo and summing outwards,i.e. ^ Msat- We 

include all post-processed subhaloes above our mass threshold of 
20 particles. As Figure|7]demonstrates at level 4 most of the finders 
(AHF, HBT, HOT6d, HSF, STF, VOBOZ) agree very well, finding 
very similar amounts of substructure both in radial location and 
mass. ROCKSTAR finds a little more structure, particularly in the 
central region where its phase-space nature works to its advantage 
and SUBFIND finds around a factor of 25 per cent less due to its 
conservative subhalo mass assignment. 

The MENDIETA finder appears to show significantly different 
results to the rest. As previously noted the ADAPTAHOP finder lo- 
cates many small subhaloes and these push up the total mass found 
in substructure above that found by the others particular in the range 
around 50-lOOkpc. We note that two of the three phase-space based 
finders (HOT6d, & HSF) have a radial performance indistinguish- 
able from real-spaced based finders. The only one to show any dif- 
ference is ROCKSTAR and it remains unclear whether or not this is 
in practice a significant improvement. 

We further like to mention (though not explicitly shown here) 
that a visual comparison akin to Figure [T|but focusing on the cen- 
tral 20 h^^kpc reveals that it is very likely that the excess mass 
found in that inner region by some of the finders such as H0T3d 
may be due to mis-identifications of the host halo as subhaloes. In 
the very central region it is difficult for the underlying real-space 
Friends-of- Friends methodology to distinguish structures from the 
background halo and so can show up either as multiple detections, 
or no structure at all. 



5 SUMMARY & CONCLUSIONS 

We have used a suite of increasing resolution models of a single 
Milky Way sized halo extracted from a self-consistent cosmo log- 
ical simulation (i.e. the Aquarius suite JSpringel et al.l 120081) ) to 
study the accuracy of substructure recovery by a wide range of 
popular substructure finders. Each participating group analysed in- 
dependently as many levels of the Aquarius-A dataset at redshift 
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Figure 6. Cumulative number count of subhaloes above the indicated i)max value within a radius of 250 kpc/h from the fiducial halo centre after standardised 
post-processing (see Section l4~ll . The arrows indicate the number of particles interior to rmax, the position of the peak of the rotation curve. The bottom plot 
shows the relative offset from the mean of the cumulative count. 



2 = as they could manage and returned lists of particles they 
associated with any subhalo they found. These lists were post- 
processed by a single uniform analysis pipeline. This pipeline em- 
ployed a standard fixed definition of the subhalo centre and subhalo 
mass, and employed a standard methodology for deriving iimax- 
This analysis was used to produce cumulative number counts of 
the subhaloes and examine how well each finder was able to locate 
substructure. 

We find a remarkable agreement between the finders which 



are based on widely different algorithms and concepts. The finders 
agree very well on the presence and location of subhaloes and quan- 
tities that depend on this or the inner part o f the halo are amazi ngly 
well and reliably recovered. We agree with Knebe et al.lfcOllh that 
^max is a good parameter by which to rank order the haloes (in 
this case subhaloes). However, we also show that as Wmax is only 
dependent upon the inner 20 per cent or less of the subhalo parti- 
cles around 100 particles are required to be within the subhalo for 
this measure to be reliably recovered. Quantities that depend on the 
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Figure 7. Cumulative plot of the enclosed mass within subhaloes as a function of radial distance from the fiducial centre of the host halo. 



outer parts of the subhaloes, such as the total mass, are still recov- 
ered with a scatter of around 10 per cent but are more dependent 
upon the exact algorithm employed both for unbinding (intrinsic 
to each finder) and to define the outer edge (given by the common 
post-processing applied here). 

The most difficult region within which to resolve substruc- 
tures is the very centre of the halo which has, by definition, a very 
high background density. In this region real-space based finders are 
expected to struggle whereas the full six-dimensional phase-space 
based finders should do better. In practice ROCKSTAR is the only 
phase-space based finder that shows any indication of this (and this 
difference becomes less pronounced as the resolution is increased); 
but we cannot rule out mis-identifications of the host halo as sub- 
haloes at this stage. We conclude that, as yet, none of the phase- 
space based finders present a significant improvement upon the best 
of the more traditional real-space based finders. 

Convergence studies indicate that identified subhaloes con- 
taining less than 100 particles tend to be under-resolved and these 
objects grow slightly in mass if a higher resolution study is used. 
This could be due to the fact that particles in the outer regions 
of these subhaloes are stripped more readily at lower resolution 



or it could be an artefact of the difficulty of measuring the po- 
tential (and hence completing any unbinding s atisfactorily') with 
this small number of particles. Several s tudies (IKase et aLluOOTI : 
IPilipenko eT^ l2009l : iTrenti et al.l 1201 Oh have indicated the un- 
reliability of halo properties (other than physical presence) for 
(sub)haloes of this size or less. 

Fitting power law slopes to the convergence studies of each 
finder indicates that the logarithmic slope of the cumulative num- 
ber count is less than 1. While this is only confirmed for a sin- 
gle halo within a single cosmology, and ignoring any mass in tidal 
streams, the result appears to be robust as it is found for all the 
high-resolution finders employed in this study. This indicates that 
the larger substructures are the most important ones and that higher 
levels of the (sub)subhalo hierarchy play a less significant dynami- 
cal role. 

We like to close with a brief note on the removal of gravi- 
tationally unbound particles for subhaloes. We have seen that the 
omission of such a procedure most certainly leads to rather dis- 
tinct results. However, we cannot convincingly deduce whether or 
not this will lead to more small mass objects (as is the case for 
ADAPTAHOP) or to objects more massive in general (also seen for 
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ADAPTAHOP), likely both will occur. But we confirm that the exact 
differences between a subhalo catalogue based upon a halo find- 
ing method with and without unbinding depend on the actual algo- 
rithm to collect the initial set of particles to be considered part of 
the subhalo: we performed an analysis of the level 4 data with AHF 
switching the unbinding part off ending up with a subhalo mass 
function that was only different at the higher mass end (not shown 
here though) as opposed to the ADAPTAHOP results; but both these 
codes differ substantially in the way of assigning the primary par- 
ticle set to a subhalo. 

It should be not ed that the Aquarius-A halo is a relatively qui- 
escent halo cWang et alj|201 ih . not having been subject to many 
mergers. Investigation of other haloes, and those produced by other 
simulation code would be interesting to compare. Therefore more 
studies focusing on the actual halo catalogues returned by each 
finder (as opposed to the particle lists used here); other cosmologi- 
cal simulations and different simulated scenarios (such as disrupted 
galaxies); the detailed analysis of sub-substructure (which is only 
really practical at level 1); and other subhalo properties such as spin 
parameter and shape, as well as more detailed resolution studies for 
those codes providing an analysis of all levels will be deferred to 
future papers. 
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